#ifndef BVPS_H_
#define BVPS_H_

#include <iostream>
#include <vector>
#include <math.h>
#include <numeric>
#include <fstream>
#include "Eigen/Dense"
#include "circle.hpp"
#include "function.hpp"
#include "norm.hpp"

using namespace std;
using namespace Eigen;

class BVPs
{
private:
    double h;
    circle c;
    Function f;
    BoundaryFunc g; // g is the boundary condition, i.e. U(x, y)
    vector<vector<point>> p_list;
    vector<double> sol;
    int m;
    MatrixXd A;
    VectorXd f_E;

public:
    /************************************************设置网格************************************************/
    void set_grid(int _n); // 网格化,第i行j列坐标为(jh, ih), i,j = 0, 1, ..., m, m+1

    /************************************************设置圆盘位置************************************************/
    void set_disk(); // 标记点是否在圆内

    /************************************************解决（不）规则边界Dirichlet问题************************************************/
    void solve_Dirichlet_bound();

    /************************************************解决规则边界Neumann问题************************************************/
    void solve_Neumann_bound();

    /************************************************输出matlab代码************************************************/
    void display();
    /************************************************计算误差************************************************/
    void compute_norm();
};

void BVPs::set_grid(int _n) // 网格化,第i行j列坐标为(jh, ih), i,j = 0, 1, ..., m, m+1
{
    m = _n;
    h = 1.0 / (m + 1);
    for (int i = 0; i <= m + 1; i++)
    {
        vector<point> vtmp;
        for (int j = 0; j <= m + 1; j++)
        {
            point tmp(j * h, i * h);
            vtmp.push_back(tmp);
            // cout << tmp.insquare;
        }
        // cout << endl;
        p_list.push_back(vtmp);
    }
    // cout << endl;
}

void BVPs::set_disk() // 标记点是否在圆内
{
    point center(0.5, 0.5);
    c.set(center, 1.0 / 16.0);
    for (int i = 0; i <= m + 1; i++)
    {
        for (int j = 0; j <= m + 1; j++)
        {
            p_list[i][j].indisk = c.is_in_disk(p_list[i][j]);
            // cout << p_list[i][j].indisk;
        }
        // cout << endl;
    }
}

void BVPs::solve_Dirichlet_bound()
{
    int dimention = pow(m + 2, 2);
    A = MatrixXd::Zero(dimention, dimention);
    f_E = VectorXd::Zero(dimention);
    int cnt = 0;
    for (int i = 0; i <= m + 1; i++)
    {
        for (int j = 0; j <= m + 1; j++)
        {
            if (p_list[i][j].insquare != 0) // 在方形边界上
            {
                f_E[cnt] = g(p_list[i][j]);
                A(cnt, cnt) = 1;
            }
            else // 在方形内
            {
                if (p_list[i][j].indisk == 1) // 还在圆形内
                {
                    f_E[cnt] = 0;
                    A(cnt, cnt) = 1;
                }
                else if (p_list[i][j].indisk == -1) // 方形内圆盘上
                {
                    f_E[cnt] = g(p_list[i][j]);
                    A(cnt, cnt) = 1;
                }
                else // 在方形内圆盘外
                {
                    double a1 = h, a2 = h, b1 = h, b2 = h;
                    f_E[cnt] = f(p_list[i][j]);
                    /******************左边是否在圆内*******************/
                    if (p_list[i][j - 1].indisk == 1)
                    {
                        point tmp = c.CrossY_RightX(i * h);
                        a1 = j * h - tmp.x;
                        f_E[cnt] += 2.0 * g(tmp) / (a1 * (a1 + a2));
                    }
                    /******************右边是否在圆内*******************/
                    else if (p_list[i][j + 1].indisk == 1)
                    {
                        point tmp = c.CrossY_LeftX(i * h);
                        a2 = tmp.x - j * h;
                        f_E[cnt] += 2.0 * g(tmp) / (a2 * (a1 + a2));
                    }

                    if (p_list[i][j - 1].indisk != 1)
                    {
                        A(cnt, cnt - 1) = -2.0 / (a1 * (a1 + a2));
                    }
                    if (p_list[i][j + 1].indisk != 1)
                    {
                        A(cnt, cnt + 1) = -2.0 / (a2 * (a1 + a2));
                    }

                    /******************上边是否在圆内*******************/
                    if (p_list[i + 1][j].indisk == 1)
                    {
                        point tmp = c.CrossX_BottomY(j * h);
                        b1 = tmp.y - i * h;
                        f_E[cnt] += 2.0 * g(tmp) / (b1 * (b1 + b2));
                    }
                    /******************下边是否在圆内*******************/
                    else if (p_list[i - 1][j].indisk == 1)
                    {
                        point tmp = c.CrossX_TopY(j * h);
                        b2 = i * h - tmp.y;
                        f_E[cnt] += 2.0 * g(tmp) / (b2 * (b1 + b2));
                    }

                    if (p_list[i + 1][j].indisk != 1)
                    {
                        A(cnt, cnt + m + 2) = -2.0 / (b1 * (b1 + b2));
                    }
                    if (p_list[i - 1][j].indisk != 1)
                    {
                        A(cnt, cnt - m - 2) = -2.0 / (b2 * (b1 + b2));
                    }
                    A(cnt, cnt) = 2.0 * (1.0 / (a1 * a2) + 1.0 / (b1 * b2));
                    // cout << i << "," << j << ":" << a1 << " " << a2 << " " << b1 << " " << b2 << " " << A(cnt, cnt) << endl;
                }
            }
            cnt++;
        }
    }
    // cout << A << endl;
    VectorXd U = A.lu().solve(f_E);
    sol.clear();
    for (int i = 0; i < dimention; i++)
    {
        sol.push_back(U(i));
    }
    display();
    compute_norm();
}

void BVPs::display()
{
    ofstream outfile;
    outfile.open("../output/test.m");
    outfile << "x = [";
    for (int i = 0; i <= m + 1; i++)
    {
        for (int j = 0; j <= m + 1; j++)
        {
            outfile << p_list[i][j].x << ",";
        }
    }
    outfile << "];" << endl;

    outfile << "y = [";
    for (int i = 0; i <= m + 1; i++)
    {
        for (int j = 0; j <= m + 1; j++)
        {
            outfile << p_list[i][j].y << ",";
        }
    }
    outfile << "];" << endl;

    outfile << "z = [";
    for (int i = 0; i < sol.size(); i++)
    {
        outfile << sol[i] << ",";
    }
    outfile << "];" << endl;
    outfile << "scatter3(x, y, z, \"*\");" << endl;
}

void BVPs::compute_norm()
    {
        int cnt = 0;
        vector<double> err;
        for (int i = 0; i <= m + 1; i++)
        {
            for (int j = 0; j <= m + 1; j++)
            {
                if (sol[cnt] != 0)
                    err.push_back(g(p_list[i][j]) - sol[cnt]);
                cnt++;
            }
        }
        Norm err_norm;
        cout << "infinite norm is:" << err_norm.inftynorm(err) << endl;
        cout << "onenorm is:" << err_norm.onenorm(err) << endl;
        cout << "twonorm is:" << err_norm.twonorm(err) << endl;
        // err_norm.output(err);
    }

#endif